%% Manylabs benchmark analysis
% Replicates Figure 2 bottom row, Figure A-2 bottom row

% Set max iteration number and tolerance
max_iter = 5000;
tol = 0.0001;

% Random seed
seed = 101;

% Sample size of simulated coverage
sim_cov = 5000;

%% Case 1: N=66, J=5 (All effects included, nu estimated by EB)
N = 66;
J = 5;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet1 = 'N=66, J=5';
estimates1 = xlsread('../Data/ManyLabs_Data.xlsx', sheet1, 'C3:H68'); 
se1 = xlsread('../Data/ManyLabs_Data.xlsx', sheet1, 'N3:S68');   

% Randomly assign baseline and validation studies
[hat_theta_1, hat_vartheta_1, rep_base_se2_1, rep_val_se2_1] = random_assignment(estimates1, se1, seed);

% Estimate nu from pooled data
EB_est_nu_1 = MLE_nu(hat_theta_1, rep_base_se2_1, 0, max_iter, tol);

% Compute and plot ecf as a function of nu
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[empirical_cov_list_1, cov_p_value_list_1] = ecf_over_nu(hat_theta_1, rep_base_se2_1, hat_vartheta_1, rep_val_se2_1, meta_nu_list, cov_stat_sim);
cov_plot_1 = plot_cov_over_nu(empirical_cov_list_1, meta_nu_list, EB_est_nu_1);
saveas(cov_plot_1, '../Results/nu_plot_N_66_J_5.png');

%% Case 2: N=132, J=2 (All effects included, nu estimated by EB)
N = 132;
J = 2;
num_blocks = 11;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet2 = 'N=132, J=2';
estimates2 = xlsread('../Data/ManyLabs_Data.xlsx', sheet2, 'C3:E134'); 
se2 = xlsread('../Data/ManyLabs_Data.xlsx', sheet2, 'K3:M134'); 

% Randomly assign baseline and validation studies
[hat_theta_2, hat_vartheta_2, rep_base_se2_2, rep_val_se2_2] = random_assignment(estimates2, se2, seed);

% Estimate nu from pooled data
EB_est_nu_2 = MLE_nu(hat_theta_2, rep_base_se2_2, 0, max_iter, tol);

% Compute and plot ecf as a function of nu
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[empirical_cov_list_2, cov_p_value_list_2] = ecf_over_nu(hat_theta_2, rep_base_se2_2, hat_vartheta_2, rep_val_se2_2, meta_nu_list, cov_stat_sim);
cov_plot_2 = plot_cov_over_nu(empirical_cov_list_2, meta_nu_list, EB_est_nu_2);
saveas(cov_plot_2, '../Results/nu_plot_N_132_J_2.png');

%% Case 3: N=48, J=5 (three irreplicable studies excluded, nu estimated by EB)
N = 48;
J = 5;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet3 = 'N=48, J=5';
estimates3 = xlsread('../Data/ManyLabs_Data.xlsx', sheet3, 'C3:H50'); 
se3 = xlsread('../Data/ManyLabs_Data.xlsx', sheet3, 'N3:S50');   

% Randomly assign baseline and validation studies
[hat_theta_3, hat_vartheta_3, rep_base_se2_3, rep_val_se2_3] = random_assignment(estimates3, se3, seed);

% Estimate nu from pooled data
EB_est_nu_3 = MLE_nu(hat_theta_3, rep_base_se2_3, 0, max_iter, tol);

% Compute and plot ecf as a function of nu
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[empirical_cov_list_3, cov_p_value_list_3] = ecf_over_nu(hat_theta_3, rep_base_se2_3, hat_vartheta_3, rep_val_se2_3, meta_nu_list, cov_stat_sim);
cov_plot_3 = plot_cov_over_nu(empirical_cov_list_3, meta_nu_list, EB_est_nu_3);
saveas(cov_plot_3, '../Results/nu_plot_N_48_J_5.png');

%% Case 4: N=96, J=2 (three irreplicable studies excluded, nu estimated by EB)
N = 96;
J = 2;
meta_nu_list = 0:0.01:0.5;

% Extract study estimates and ses
sheet4 = 'N=96, J=2';
estimates4 = xlsread('../Data/ManyLabs_Data.xlsx', sheet4, 'C3:E98'); 
se4 = xlsread('../Data/ManyLabs_Data.xlsx', sheet4, 'K3:M98');   

% Randomly assign baseline and validation studies
[hat_theta_4, hat_vartheta_4, rep_base_se2_4, rep_val_se2_4] = random_assignment(estimates4, se4, seed);

% Estimate nu from pooled data
EB_est_nu_4 = MLE_nu(hat_theta_4, rep_base_se2_4, 0, max_iter, tol);

% Compute and plot ecf as a function of nu
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[empirical_cov_list_4, cov_p_value_list_4] = ecf_over_nu(hat_theta_4, rep_base_se2_4, hat_vartheta_4, rep_val_se2_4, meta_nu_list, cov_stat_sim);
cov_plot_4 = plot_cov_over_nu(empirical_cov_list_4, meta_nu_list, EB_est_nu_4);
saveas(cov_plot_4, '../Results/nu_plot_N_96_J_2.png');
